# Authors: Emil & Ivetta

# this script does main conjoint analysis
# beware cjplot package needs categorical variables to be factors, otherwise it doesn't understand
# this script produces plots: 
# Figure 2 in the main text
# Tables 8, 9 and Figure 5 in Appendix H
# Figures 6 and 7 in Appendix I
# Table 10 and 11 in Appendix J

pacman::p_load(haven, dplyr, labelled, xlsx, cjoint, sandwich, lmtest, cregg,
               ggtext, glue, ggsci)
#install.packages("ggstance")
library(cregg)
#detach("package:rlang", unload = TRUE)


# For afcp estimand
#remotes::install_github("astrezhnev/afcp", build_vignettes = TRUE)
library(afcp)

#setting cool ggtheme
theme_Publication <- function(base_size=18, base_family="Helvetica") {
  library(grid)
  library(ggthemes)
  (theme_foundation(base_size=base_size, base_family=base_family)
    + theme(plot.title = element_text(face = "bold",
                                      size = rel(0.9), hjust = 0.5),
            text = element_text(),
            panel.background = element_rect(colour = NA),
            plot.background = element_rect(colour = NA),
            panel.border = element_rect(colour = NA),
            axis.title = element_text(face = "bold",size = rel(0.8)),
            axis.title.y = element_text(angle=90,vjust =2),
            axis.title.x = element_text(vjust = -0.2),
            axis.text = element_text(), 
            axis.line = element_line(colour="black"),
            axis.ticks = element_line(),
            panel.grid.major = element_line(colour="#f0f0f0"),
            panel.grid.minor = element_blank(),
            legend.key = element_rect(colour = NA),
            legend.position = "bottom",
            legend.direction = "horizontal",
            legend.key.size= unit(0.2, "cm"),
            legend.margin = unit(0, "cm"),
            legend.title = element_text(face="italic"),
            plot.margin=unit(c(10,5,5,5),"mm"),
            strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
            strip.text = element_text(face="bold")
    ))
  
}



#import conjoint df enriched by all the vars from the second wave survey
df <- read_sav("drafts/conjoint_solidarity/perspectives_docs/revision_2/replication_code/input_data.sav")
df <- as_factor(df)




#data frame that contains only those who are beyond Russia
df_mig <- df %>%
  filter(location %in% "Да") %>%
  filter(when_left %in% "Да") #filter out old emigrants



# to make it a new line on the plot we need to change \n onto <br />

library(forcats)

# Assuming df_mig is your dataframe and Motivation is the factor variable
df_mig <- df_mig %>%
  mutate(Motivation = fct_recode(Motivation,
                                 `Can't accept Russian politics +<br />Was arrested in rallies` = "Can't accept Russian politics +\nWas arrested in rallies"))


# AMCE


amce <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Ethnicity + Motivation, 
         id = ~ResponseId, 
         estimate = "amce")




# Creating caps lock for attributes in plots
label_rename <- c(
  Age = "Age:",
  Gender = "Gender:",
  Children = "Children:",
  Profession = "Profession:",
  Ethnicity = "Ethnic background:",
  Motivation = "Motivation for emigration:"
)

# Assuming the dataframe is stored directly in the amce object (adjust accordingly if it's nested further)
amce$feature <- ifelse(amce$feature %in% names(label_rename), 
                       label_rename[amce$feature], 
                       amce$feature)


# Define labels to be bold
bold_labels <- unique(amce$feature)
bold_labels <- paste0(bold_labels, collapse = "|")

# Custom function to apply bold to specific labels
highlight <- function(x, pat, color="black", family="") {
  # Remove parentheses
  x <- gsub("[()]", "", x)
  
  # Apply bold to labels that match the pattern
  ifelse(grepl(pat, x), glue::glue("<b style='font-family:{family}; color:{color}'>{x}</b>"), x)
}


# Full version for the appendix
amce_plot <- plot(amce, vline = 0, size = 2, legend_pos = "none")+ 
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none") + 
  geom_point(size = 2) +
  scale_color_nejm() +
  scale_y_discrete(labels= function(x) highlight(x, bold_labels, "black")) +
  theme(axis.text.y=element_markdown())





ggsave("drafts/conjoint_solidarity/plots/main_plots/amce.png",
       amce_plot,
       dpi = 500,
       width = 5, height = 4.5)

# Shorter version (only relevant attributes) for the main text
amce <- amce %>% 
  filter(!feature %in% c("Gender:", "Age:", "Children:", "Profession:"))

amce_plot <- plot(amce, vline = 0, size = 2, legend_pos = "none")+ 
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none") + 
  geom_point(size = 2) +
  scale_color_nejm() +
  scale_y_discrete(labels= function(x) highlight(x, bold_labels, "black")) +
  theme(axis.text.y=element_markdown())



ggsave("drafts/conjoint_solidarity/plots/main_plots/amce_short.png",
       amce_plot,
       dpi = 500,
       width = 5, height = 4)


# Add multiple comparsion corrections

amce <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Ethnicity + Motivation, 
           id = ~ResponseId, 
           estimate = "amce")

library(stats)

options(scipen = 999)

#create a named vector of p-values from cregg output
p_named <- na.omit(amce$p)
names(p_named) <- amce$level[!is.na(amce$p)]

p_corrected_bh <- p.adjust(p_named,
         method = "BH") %>% round(5) # all effects holds

p_corrected_bh <- data.frame(level = names(p_corrected_bh), 
                             'p_fdr' = p_corrected_bh, 
                             stringsAsFactors = FALSE,
                             row.names = NULL)

p_corrected_bonf <- p.adjust(p_named,
         method = "bonferroni") %>% round(5)

p_corrected_bonf <- data.frame(level = names(p_corrected_bonf), 
                             'p_bonf' = p_corrected_bonf, 
                             stringsAsFactors = FALSE,
                             row.names = NULL)

# assign these values to amce data frame along levels

amce <- left_join(amce,
                  p_corrected_bh)
amce <- left_join(amce,
                  p_corrected_bonf)


options(scipen = 0)

# Remove unecessary columns for table
amce <- amce %>% 
  select(-outcome, -statistic, -lower, -upper) %>%
  #renaming columns
  rename(Feature = feature,
         Level = level,
         Estimate = estimate,
         "Std. error" = std.error,
         "z-score" = z,
         "p-value" = p,
         "FDR corr." = p_fdr,
         "Bonferroni\ncorr." = p_bonf)




table_output <- xtable::xtable(amce, 
                               caption = "\\label{tab:amce_main_effects}Average Marginal Component Effects of potential emigrants attributes' on the probability of receiving Help. FDR (p_fdr) and Bonferroni (p_bonf) corrected p-values added.", 
                               include.rownames = FALSE,
                               digits = 3)

print(table_output, include.rownames=FALSE, include.colnames=TRUE)


# Marginal Means

mm <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Ethnicity + Motivation, 
         id = ~ResponseId, 
         estimate = "mm")
# Creating caps lock for attributes in plots
label_rename <- c(
  Age = "Age:",
  Gender = "Gender:",
  Children = "Children:",
  Profession = "Profession:",
  Ethnicity = "Ethnic background:",
  Motivation = "Motivation for emigration:"
)

mm$feature <- ifelse(mm$feature %in% names(label_rename), 
                       label_rename[mm$feature], 
                     mm$feature)


mm_plot <- plot(mm, vline = 0.5, size = 2, legend_pos = "none") + 
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none") +
  # theme(
  #   plot.background = element_rect(color = "black", fill = NA, size = 1.5, linetype = "solid"), # Add black border around the plot
  #   panel.background = element_blank() # Ensure that the panel background doesn't overlap with the plot border
  # ) +
  geom_point(size = 2) +
  scale_color_nejm() +
  scale_y_discrete(labels= function(x) highlight(x, bold_labels, "black")) +
  theme(axis.text.y=element_markdown())

ggsave("drafts/conjoint_solidarity/plots/main_plots/mm.jpeg",
       mm_plot,
       dpi = 500,
       width = 5, height = 4.5)

mm <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Ethnicity + Motivation, 
         id = ~ResponseId, 
         estimate = "mm")

mm <- mm %>% 
  select(-outcome, -statistic, -lower, -upper) %>%
  #renaming columns
  rename(Feature = feature,
         Level = level,
         Estimate = estimate,
         "Std. error" = std.error,
         "z-score" = z,
         "p-value" = p)


table_output <- xtable::xtable(mm, 
                               caption = "\\label{tab:mm_main_effects}Marginal means of potential emigrants attributes'.", 
                               include.rownames = FALSE,
                               digits = 3)

print(table_output, include.rownames=FALSE, include.colnames=TRUE)


############################################################
# average feature choice probability (AFCP)
# for afcp package we need to recode profile variable, now it has values 1.1,1.2,2.1,2.2, and the function takes it as four profiles instead of just two for different tasks
############################################################

df_mig <- df_mig %>%
  mutate(profile_2 = case_when(profile %in% c("1.1", "2.1") ~ 1,
                               profile %in% c("1.2", "2.2") ~ 2))

levels(df_mig$Motivation) # function afcp() seems to dislike <br />, we need to remove it

df_mig <- df_mig %>%
  mutate(Motivation = fct_recode(Motivation,
                                 "Was arrested in rallies" = 
                                   "Can't accept Russian politics +<br />Was arrested in rallies"
  ))


table(df_mig$Motivation)
table(df_mig$Ethnicity)


levels(df_mig$Motivation)
levels(df_mig$Ethnicity)



library(cjoint)
# Motivation
amce_results <- cjoint::amce(choice ~ `Motivation`,
                             data = df_mig, cluster = TRUE, respondent.id = "ResponseId", design = "uniform")


afcp_motiv <- afcp(amce_results, 
             respondent.id = "ResponseId", 
             task.id = "task",
             profile.id = "profile_2",
             attribute = "Motivation")


afcp_motiv$afcp
afcp_motiv$wald

# Ethnicity
amce_results <- cjoint::amce(choice ~ `Ethnicity`,
                     data = df_mig, cluster = TRUE, respondent.id = "ResponseId", design = "uniform")


afcp_ethnic <- afcp(amce_results, 
           respondent.id = "ResponseId", 
           task.id = "task",
           profile.id = "profile_2",
           attribute = "Ethnicity")

afcp_ethnic$afcp
afcp_ethnic$wald


# Profession
amce_results <- cjoint::amce(choice ~ `Profession`,
                             data = df_mig, cluster = TRUE, respondent.id = "ResponseId", design = "uniform")


afcp_proff <- afcp(amce_results, 
             respondent.id = "ResponseId", 
             task.id = "task",
             profile.id = "profile_2",
             attribute = "Profession")

afcp_proff$afcp
afcp_proff$wald



# Now making it one:

# Add category and bind AFCP tables
afcp_all <- bind_rows(
  afcp_motiv$afcp  %>% mutate(attribute = "motive"),
  afcp_proff$afcp  %>% mutate(attribute = "profession"),
  afcp_ethnic$afcp %>% mutate(attribute = "ethnic")
) %>%
  relocate(attribute) %>%
  select(-conf_high, -conf_level, -conf_low)

# Add category and bind Wald tables
wald_all <- bind_rows(
  afcp_motiv$wald  %>% mutate(attribute = "motive"),
  afcp_proff$wald  %>% mutate(attribute = "profession"),
  afcp_ethnic$wald %>% mutate(attribute = "ethnic")
) %>%
  relocate(attribute)

library(knitr)

kable(afcp_all, format = "latex", booktabs = TRUE, digits = 3, caption = "AFCP Results with Category")
kable(wald_all, format = "latex", booktabs = TRUE, digits = 3, caption = "Wald Test Results with Category")



####################################################
##############Interactions (ACIE)###################
####################################################


# Interaction with profession
acie_4 <- cj(df_mig, choice ~ Age + Gender + Children + Ethnicity + Motivation, id = ~ResponseId, estimate = "amce", 
             by = ~Profession)

diff_acie_4 <- cj(df_mig, choice ~ Age + Gender + Children + Ethnicity + Motivation, id = ~ResponseId, estimate = "amce_diff", 
                  by = ~Profession)

acie_plot_profession <- plot(rbind(acie_4, diff_acie_4)) + ggplot2::facet_wrap(~BY, ncol = 3L)


ggsave("drafts/conjoint_solidarity/plots/main_plots/acie_profession.png",
       acie_plot_profession,
       dpi = 500,
       width = 8, height = 6)

# Interaction with Ethnicity
acie_5 <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Motivation, id = ~ResponseId, estimate = "amce", 
             by = ~Ethnicity)

diff_acie_5 <- cj(df_mig, choice ~ Age + Gender + Children + Profession + Motivation, id = ~ResponseId, estimate = "amce_diff", 
                  by = ~Ethnicity)

acie_plot_ethnic <- plot(rbind(acie_5, diff_acie_5)) + ggplot2::facet_wrap(~BY, ncol = 3L)

ggsave("drafts/conjoint_solidarity/plots/main_plots/acie_ethnic.png",
       acie_plot_ethnic,
       dpi = 500,
       width = 8, height = 6)




